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Abstract 

For the contact of two finite portions of interacting rigid crystalline surfaces, 
we compute the the pinning energy barrier dependency on the misfit angle 
and contact area. This simple model allows us to investigate a broad contact- 
size and angular range, thus obtaining the statistical properties of the energy 
barriers opposing sliding for a single asperity. These data are used to generate 
the distribution of static frictional thresholds for the contact of polycrystals, as 
in dry or even lubricated friction. This distribution is used as the input of a 
master equation to predict the sliding properties of macroscopic contacts. 

Keywords: tribology, mechanical contacts, atomic scale friction, boundary 
lubrication, relative crystalline orientation, size scaling 



1. Introduction 

In a regime of dry friction or boundary lubrication of a single contact, such 
as studied by atomic force microscope (AFM) techniques, the friction force 
depends substantially and nontrivially on the relative crystalline orientation of 
the facing surfaces, as was demonstrated experimentally by Dienwiebel et al. 
Special angles lead to superlubric sliding, but tend to be energetically unfavor- 
able Q. Depending on the contact mechanical details and the sliding speed, 
such superlubric orientations could have long enough time to reconstruct, ap- 
proaching a state with a lower free energy but characterized by a higher barrier 
(aging) , or be retained for long enough for them to provide a substantial super- 
lubric contribution to the overall sliding dynamics. The connection between the 
nanoscale, where friction occurs through the breaking and formation of local 
contacts, and the meso- or macroscale, where many breaking junctions interact 
elastically, is commonly described by ear thq uake-type models 
or by a master equation approach [Io[ HH, 12 1, or by models inspired to 



the Greenwood- Williamson one [13, [lj] such as the sub-boundary lubrication 
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model 15|, [lg, [l2 1- Except the case of ideally flat surfaces such as mica in 



surface- force apparatus (SFA) experiments, contact is always realized at micro- 
scopic asperities, whose size ranges typically in the nanometer to micrometer 
range. Even when a lubricant is present at the contacting asperities, in a bound- 
ary or sub-boundary lubrication context, the residual lubricant is often frozen 
by pressure and confinement, and hydrodynamic viscous sliding is replaced by 
a static (stick-slip) contact-breaking regime, which we focus on in the present 
work. In the multiasperity contact, where relative orientation of individual 
asperities is not really under control, the most important information to be ex- 
tracted from a single- asperity model is the probability distribution P c {e a ) of the 
slip activation barriers e a . 

Molecular dynamics (MP) approaches to lubricated sliding friction [l8l [Til 
[H ED, HI, 0, 0, [H, HI [23, Q are usually forced to use some form of pe- 
riodic boundary conditions (PBC) in order to prevent the escape of lubricant 
atoms (molecules) from the contacting region under high load, and to keep the 
simulation size under control. PBC mimic the infinite size limit, which might 
not be especially appropriate for sub-micrometer-size contacts between sliding 
surfaces. Moreover, a fully atomistic model would be computationally too de- 
manding for a full statistical study of the size and angular dependency of the 
characteristics of contacts. 

To study the contact, we introduce a simple rigid model for a finite-size 
breaking junction realized by the contact of a finite portion of two different 
crystalline surfaces. Such a rigid model could not possibly account for wear 
or for the dissipation occurring at contact breaking, as could MD simulations 
instead, but provides semi-quantitative estimates of the barriers opposing the 
onset of sliding, i.e., the static friction thresholds. The simple model allows 
us to evaluate the relevant statistical distribution of barrier energies. This 
distribution enters as a basic ingredient in the master-equation formulation, 
which, depending on the general shape of this distribution leads to different 
macroscopic sliding regimes, either stick-slip or smooth-sliding dynamics [plliol]. 

The present paper reports progress beyond our MD study of Ref . [28] . For 
the MD simulations of that work, we used a fixed size of the contact area, and 
applied PBC to impose a given commensurability ratio with minimal boundary 
effects, and to obtain a fair comparison of different misfit angles (f>. Of course, 
the friction force /, and in particular its dependency on the angle 4>, does change 
with the system size. For a larger contact area, increasingly many incommen- 
surate angles emerge, producing a more irregular dependency /(</>). In the limit 
of infinite size, f(4>) develops an infinite number of singularities [29(, but this 
limit is not to be taken too seriously, since in practice a single contact of a given 
crystalline orientation never exceeds a fraction of /1m 2 in practical experimental 
conditions. 

A numerical determination of the size dependency of f{4>) would constitute a 
formidable task for MD simulations, since larger sizes are not only individually 
more expensive computationally, but would require a finer and finer sampling of 
angles and longer equilibration times. It is instead straightforward to obtain 
the size dependency of the rigid-island model. Thus, the main goal of the present 



2 



work is to find, at least qualitatively, the shape of the distribution of static 
thresholds P c (e a ), which is the main input of the master- equation approach to 
friction on meso/macro-scale, for a contact of polycrystalline surfaces. 

In Sec. [3] we spell out the details of this rigid-contact model for the analy- 
sis of the static friction barrier realized by the contact of a crystalline surface 
with the boundary lubricant layer, which we assume in a close-packed ordered 
configuration. Sections [3] and [4] discuss the angle and size dependency of the 
friction of a nanocontact, and we compare our results with those of the MD 
lubricated model |28[. The basic implications for macroscopic sliding are dis- 
cussed in Sec. [5] within a master-equation approach. Section [5] summarizes and 
discusses our results. 



2. The rigid-island model 

We represent the sliding contact, or the solidified boundary lubricant at the 
contact, with a finite rigid crystalline layer of lattice constant a, consisting of 
N point-like atoms. We put it in interaction with a substrate potential which 
is also rigid and periodic, e.g., a sinusoid of a generally different period a s and 
amplitude Vq. In the case of a one-dimensional system, one can easily find an 
explicit expression for the activation energy barrier for the onset of motion along 
the chain: 

sin(27r Na/a s 



V 



(1) 



sin(27ra/a s ) 

Accordingly, for suitable values of the lattice constant, namely for a = n a s /(2N) 
with integer n (but not a multiple of N) , this activation barrier vanishes and the 
chain moves freely. For a nonrigid layer, the activation energy remains nonzero 
for all values of a, but still reaches the first minimum at a/a s oc 1/N: the motion 
in such a case is of a "caterpillar" type (for details see Refs. 3Ct l31|). 

A two-dimensional "lubricant" island advancing over the 2D periodic sub- 
strate should exhibit a similar behavior: in particular, a minimum for the acti- 
vation energy is expected at a/a s oc 1/y/N. The 2D system, however, has one 
extra degree of freedom, the rotation. We expect that the activation energy 
would achieve minima for specific misfit angles. 

To investigate this pattern of minima, we consider a rigid island of size ./V 
with a triangular lattice interacting with the simplest 2D substrate periodic 
potential 

V(x,y) = V (sins + sin y) (2) 

of square symmetry and lattice spacing a s — 2it. The atomic coordinates of the 
approximately square-shaped island are Xij = X + (i+j/2)a and yi j = Y + jh, 
where h = ay/3/2, and the indexes % — —j/2, . . . , —j/2+m—l, j = 0, . . . , rij — 1. 
The number of atoms in the rigid flake is N — riiUj, with n^a njh. X 
and Y are the degrees of freedom for the translation of the rigid island. If 
we rotate the island by an angle </>, then the atomic coordinates change to 
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x i,j = cos 4> — Vi,j sin </> and yij — Xij sin (j> + jjij cos <j). For a fixed misfit 
angle </>, the total potential energy of the island in contact with the substrate is 

U(X,Y)=J2 E V ^hVi,j)- (3) 

j=0 t=— i/2 

The set of stationary points of U(X, Y), defined by 

dU/dX = and dU/dY = , (4) 

consist of four elements: one minimum U m , two saddle points U a \ and U S 2, and 
one maximum. The nature of the stationary points can be verified by computing 
the Hessian of U(X, Y) at each stationary point, but it is in fact simply provided 
by the values of the function at the points. The activation energy for sliding is 

£a = U S - U m , (5) 

where the lower saddle energy U s = min(Z7 s i, U S 2) is considered. The calcula- 
tion of the energy barrier e a is extremely fast and efficient, since it only requires 
the search of the stationary points in a 2D space, by solving a simple numerical 
equation ^ . Given the simplicity of this model, it allows us to compute e a for a 
very fine sampling of angles and contact sizes, as would be practically impossible 
if a fully atomistic simulation is used. Note that the computed barrier height is 
relevant irrespective of the direction in which the rigid island is pushed forward. 

Although not completely equivalent, the activation energy barrier is indeed 
related to the threshold force f s necessary to initiate sliding. To compare the 
results of the present rigid-island approach with the fully atomistic simulations 



of Ref. [28( | . we take a/a s = 4.14/3, and assume that f s = Ke a /a s , where k is 
a factor of order unity. This comparison is shown in Fig. [T]for an intermediate 
island size. We see that, at least qualitatively, MD and the present simple rigid 
model agree on predicting static-friction minima near a set of "special" angles. 
The singularities at the optimal angles are exact zeroes for the rigid model: here 
the two lowest stationary points of U(X, Y) mutate into a degenerate trough. 
Unfortunately, within the PBC setup of MD it is impossible to cover such a fine 
sampling of misfit angles as allowed by the rigid model one: it is thus impos- 
sible to verify to what extent the two models agree or disagree on the specific 
superlubric angular orientations. The purpose of this comparison is purely to 
highlight a general qualitative similarity of the outcomes of two models, which 
are distinct and address significantly different physics. The main difference be- 
tween the rigid model and the deformable lubricant film is that the deformable 
model has a finite barrier against sliding for all angles, even at optimal ones. 



3. Barrier versus misfit angle 

Figure [2] shows the angular dependency of the barrier against sliding £ a (<^>), 
for different sizes of the sliding island. The barrier reaches its first minimum at 
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Figure 1: (Color online) (a) The static friction force (immediately before slip) as a function 
of the misfit angle <j> for a one-layer lubricant film at zero temperature for different driving 
speeds, computed with MD simulations with in-plane PBC 1281 . The friction force and driving 
speed are given in natural model units, of the order 1 nN and ~ 1 m/s respectively. For full 
details on these simulations, see Ref. [28| |. (b) The activation energy e a (units of Vb) for the 
sliding of a rigid island composed of N = 81 = 9 X 9 lubricant atoms over the rigid square 
substrate as a function of <f>. 
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misfit angle, (/) (deg.) 

Figure 2: The activation energy e a (units of Vq) as a function of the misfit angle <j> for three 
rigid lubricant islands composed by N = 9, 156 and 896 atoms. 

a misfit angle </> m ~ 7r/ (4yN), which moves to smaller and smaller angle as the 
size N of the island grows. 

Figure [3] shows that high barrier energies e a correlate well with rather stable 
configurations characterized by a low minimum energy U m , while low-barrier su- 
per lubric angles are usually characterized by unstable (high U m ) configurations. 
A higher stability of a given angular configuration could make that configura- 
tion more likely, if the asperity is free to rotate. We take this correlation into 
account when we evaluate the distribution of depinning thresholds in the rigid 
model, based on the evaluation of e a over a fine grid of angles, as illustrated 
by Fig. |4j For a given contact size N, the distribution of activation barriers 
exhibits weak divergences produced by the round maxima of £ a (<p), plus jump 
discontinuities produced by the "kinky" maxima associated to a crossing of the 
saddle points, as illustrated for two sizes by the comparison of Figs. [3] and [4] 

If the individual contacts are allowed the freedom and a sufficiently long 
time to rotate, thermal fluctuations will lead to geometric relaxation, eventu- 
ally leading to an appropriate angular distribution P<p((t>); if one can neglect the 
interaction of the contacting grain with the rest of the slider, this distribution 
should match a Boltzmann distribution P^{4>) oc exp[— U m ((j))/kBT] of the fully 
relaxed energy U m of the grain-substrate interaction. If, on the contrary, mis- 
fit angles are frozen by the microcrystalline nature of the surfaces in contact 
for much longer than the time of the experiment, all angles are equally likely 
and P<f>(4>) is a constant (equivalent to the limiting Boltzmann distribution for 
T — > oo). Averaging with these two different weight patterns leads to related but 
significantly different distributions, as illustrated by the comparison of dashed 
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Figure 3: (Color online) A direct comparison of the energy U m of the minimum (dashed) with 
the activation barrier e a (solid), as functions of the misfit angle <f> for two sizes of the rigid 
island. Energies are in units of the substrate potential corrugation Vq. 
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Figure 4: (Color online) The distribution of activation barriers e a for two island sizes reported 
in the corresponding panels of Fig. [3] The solid line is computed assuming all misfit angles 
4> are equally likely, the dashed line weights different angles according to a Boltzmann distri- 
bution of the corresponding minimum energy U m , for a temperature kgT = Vb matching the 
typical potential barrier of a single lubricant particle. 
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Figure 5: (Color online) The probability distribution of static thresholds for the rigid domains 
averaged over the domain size N with the Poissonian distribution of average domain size (a) 
N = 80 and (b) N = 20. Inset (c): the probability distribution computed by averaging over 
an exponential size distribution exp(— N/N) with N = 50. Averaging over the misfit angle 
(j> is carried out for solid (all angles are presented equally likely) and dashed (different angles 
are weighted according to the Boltzmann distribution) lines as in Fig. f4] 

and solid lines in Fig. [4] Observe in particular that the effect of the Boltzmann 
weights is to suppress the probability of small activation barriers e a - this is a 
consequence of the stable angles (minima of U m ) being typically associated to 
high barriers e a , as remarked above (see Fig. [3]). If the atomic layer represents a 
frozen lubricant, then one should beware of other ^-dependent energy contribu- 
tions to be added to U m due to the interaction with the crystalline anisotropy 
of the asperity region of the upper slider. These extra terms would of course 
influence the Boltzmann weights in the fast-rotating condition, in a way which 
could only be predicted in a condition where the details of this interaction and 
relative crystalline alignment were given. 

4. The distribution of static thresholds 

To describe friction in a meso- or macroscopic multi-contact regime it makes 
sense to assume a distribution of contact sizes N, and obtain the statistical 
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% =£ J <£ a > 

Figure 6: (Color online) The probability distribution of static thresholds Pc(x) as a function 
of renormalized barrier heights \ = £a/(£a) for the rigid domains averaged over the domain 
size N with an exponential size distribution exp(— iV/iV) with N = 20 (dotted/black), 50 
(dot-dashed/rcd) and 80 (dashed/blue). 

contact properties by averaging over N. For the size distribution one may take 
a Poisson distribution P(N) = (N N /N\) exp(— N) with an average domain size 
N. We include clusters of nearly squared shape only. 

By combining this size distribution with the distributions of unpinning bar- 
riers for individual sizes N, calculated in the previous Section, we obtain a 
global distribution of barrier heights. The resulting distribution is displayed in 
Fig. [5l where we have smoothened the singularities of the distribution for each 
individual size N by means of the convolution with a Gaussian of full width at 
half maximum matching the average inter-peak spacing (varying from 0.1 Vq to 
2.5 Vq). We see that this distribution decays roughly exponentially by approx- 
imately two decades, and then drops rapidly due to the fast large-size decay 
of the Poisson function combined with the decreasing probability of barriers of 
increasing height, as illustrated in Fig. HJd. The choice of a Poisson distribution 
is not especially critical: a similar distribution of static thresholds is obtained 
if we assume an exponential size distribution P(N) — N^ 1 exp(— N/N), see 
Fig. [5b. 

The comparison of panels (a) and (b) of Fig.[5]shows that the average contact 
size N affects the quantitative detail of the distribution, but not its qualitative 
shape. Moreover, if we plot the individual distributions as functions of the 
dimensionless rescaled activation energy barrier 

X = £a/(e a ) , (6) 

all distributions may be roughly approximated by the exponential function 
Pc{x) — ex P(~ x) as demonstrated in Fig. [6] Remarkably, the general shape 
of the distribution of barriers of the rigid-island model resembles the distribu- 
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Figure 7: (Color online) The friction force F(X) (normalized on the X — ¥ oo limiting value 
-Ffc, corresponding to the kinetic friction force in the smooth-sliding regime) as a function of 
the displacement X of the rigid slider for different distributions of static thresholds shown in 
inset. Inset: the probability distribution of static thresholds for the rigid domains averaged 
over the domain size N with an exponential size distribution exp(— N/N) with N = 50, when 
averaging over the misfit angle <j> is weighted according to a Boltzmann distribution of the 
domain energy, for different temperatures. 



tion, shown in Fig. 9 of Ref. 28], of static thresholds in the lubricated model 
based on PBC and a single size N = 80. The main visible difference is that, 
for the deformable domains of Ref. 28], the probability of small barriers is 
significantly lower than for the rigid slider at hand. 

Small barriers s a become also suppressed when the islands can rotate so that 
angles are distributed thermally. As temperature decreases, the distribution of 
static thresholds exhibits more and more pronounced local maxima at values 
corresponding to minima of the domain's potential energy (see inset of Fig.[7])- 



5. Consequence for macroscopic sliding 

Once the distribution of static thresholds is known, we can predict the dy- 
namics of the tribological system with the help of the master-equation approach 
based on an earthquake-like model [t| [HI EH • Consider the contact of two rough 
unlubricated substrates (top and bottom) on a meso- or macroscale. If the posi- 
tion of the bottom substrate is kept fixed and the top substrate is displaced by a 
distance X, the interface force (shear stress) begins to grow, F cx X. Then, the 
domains where the stress exceeds a corresponding threshold value, start to slide, 
thus relaxing the local stress, and the increase of the total force F will degrade. 
The overall average dependence F(X) follows from a solution of this master 
equation, where the distribution of static thresholds is the input parameter, as 
discussed in detail in Refs. P. [Tol [ill] . 



11 



For the threshold distributions calculated above, typical dependences F(X) 
are shown in Fig. [JJ When the misfit angles are distributed equally likely so that 
the distribution of static thresholds is nonzero down to zero threshold and P c (e a ) 
has no sharp maxima, the force F increases monotonically with X, approaching 
the kinetic- friction force Fk characteristic of smooth sliding for X — >• oo. Thus, 
in such kind of macroscopic slider, static and kinetic friction forces coincide, and 
the motion always corresponds to smooth sliding. In contrast, when the thresh- 
old distribution exhibits well pronounced sharp maxima, like for thermalized 
domains or for an ordered homogeneous thin lubricant film, the function F(X) 
reaches a maximum (which represents the macroscopic static friction force F s ) 
greater than Fk, and then decreases as the asperities give way collectively. For 
a soft enough slider, the inequality dF(X)/dX < for some X leads to the 
appearance of the elastic instability in the system dynamics 0, [l(J ■ Under such 
conditions, the macroscopic slider may exhibit stick-slip motion, provided the 
slider is soft and the delay in contact reformation is taken into account 
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6. Discussion and Conclusion 

The calculations within the simple model at hand provide, first of all, rel- 
evant insight for a single-asperity microscopic system, e.g. for a flake sliding 
over an atomically flat surface, where the slider may rotate and search for a 
local minimum of the potential energy. As a consequence, even if sliding starts 
off in a low- friction state (e.g., in the superlubric state associated to a highly 
incommensurate <j> the flake will eventually rotate and spend most of its 
time near a local minimum (Fig. [3] indicates that for the triangle-on-square ge- 
ometry of the present model such minimum need not be <f> = 0). Accordingly, 
after a relaxation time typical of the flake rotation, friction should increase (as 
predicted by the low-£ a side drop of the distribution of Fig. [SJ. Such a behavior 
was observed experimentally and in MD simulation Observe also that an 
increased rate of thermally activated jumps across the pinning barriers would 



additionally lead to a thermolubric regime [32|, |33 1 . 

More than single-asperity experiments, the focus of the present work con- 
cerns meso- and macroscopic sliding friction. At the nanoscopic level, the fric- 
tion force produced by a sliding contact depends substantially and nontrivially 
on the relative crystalline orientation of the facing surfaces. In the present work 
we provide a basic tool to connect between the nanoscale, where friction occurs 
through the breaking and formation of local contacts, and the meso/macro-scale, 
where many breaking junctions interact elastically, as commonly described by 
an earthquake-type model or by a master-equation approach. The quantity 
that summarizes the information obtained by averaging over all possible con- 
tact sizes and angles is a probability distribution P c {e a ) of the slip activation 
barriers e a . Our simple model permits us to evaluate such a distribution of bar- 
rier energies, reaching beyond the small sizes and few rotation angles allowed by 
detailed microscopical MD simulations. This distribution is a basic ingredient 
for the master-equation formulation, which, depending on the actual shape of 
this distribution can lead to different general macroscopic sliding regimes. The 
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analysis of the shape of this distribution allows one to understand the physics of 
the meso-macroscopic sliding in terms of the underlying microscopic junction- 
breaking statistical properties. 

Two basic regimes of macroscopic sliding emerge from this model: (i) When 
supcrlubric alignments are suppressed by aging to thermodynamically more 
favorable alignments, a nonmonotonic peaked distribution P c {e a ) of barrier 
heights is obtained, which tends to induce a macroscopic stick-slip regime, (ii) In 
contrast, when the probability of weak activation barriers is sufficiently large to 
produce a monotonically decaying distribution P c (s a ), then macroscopic smooth 
sliding is possible, even in the presence of microscopic breaking-junction dynam- 
ics. 

The present simple and very idealized model is not meant to address any 
specific properties of a well-defined contacting system, but it focuses on the 
possibility to extract macroscopic statistical information out of the mechanical 
properties of contacts. Many details of real contacts are left out, including 
surface curvature, wear, local thermal expansion. For this reason, it would 
be interesting (although extremely expensive numerically) to attempt a similar 
statistical method using the MD simulations of a specific contact described in 
terms of realistic force fields and curved surfaces. While the quantitative detail 
of P c (sa) is likely to depend on the specific contacting materials, its general 
properties should mostly follow those determined by means of the present simple 
model. 
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